Method for characterising the anisotropy of the texture of a digital image

ABSTRACT

This characterizing method comprises:
         estimating ( 28 ) the scalar coefficients τ m  of an even function τ(θ) defined in [0; 2π] that minimizes the following criterion C:       

     
       
         
           
             C 
             = 
             
               
                 ∑ 
                 
                   j 
                   = 
                   1 
                 
                 
                   n 
                   i 
                 
               
                
               
                 
                   ( 
                   
                     
                       β 
                       j 
                     
                     - 
                     
                       τ 
                       * 
                       
                         Γ 
                          
                         
                           ( 
                           
                             α 
                             j 
                           
                           ) 
                         
                       
                     
                   
                   ) 
                 
                 2 
               
             
           
         
       
     
     where:
         β j  are terms estimated from an acquired digital image,   τ(θ) is a π-periodic function defined in the interval [0; 2π],   Γ(θ) is the function defined by the following relationship:       

     
       
         
           
             
               Γ 
                
               
                 ( 
                 θ 
                 ) 
               
             
             = 
             
               
                 ∫ 
                 
                   ℝ 
                   + 
                 
               
                
               
                 
                   
                      
                     
                       
                         v 
                         ^ 
                       
                        
                       
                         ( 
                         
                           ρ 
                            
                           
                               
                           
                            
                           θ 
                         
                         ) 
                       
                     
                      
                   
                   2 
                 
                  
                 
                   ρ 
                   
                     
                       
                         - 
                         2 
                       
                        
                       H 
                     
                     - 
                     1 
                   
                 
                  
                 d 
                  
                 
                     
                 
                  
                 ρ 
               
             
           
         
       
         
         
           
             
               
                 where:
               {circumflex over (V)} is the discrete Fourier transform of a convolution kernel v,   H is an estimated Hurst exponent of the acquired image,
 
f) then, calculating ( 30 ), depending on the estimate of the scalar coefficients τ m , an anisotropy index that characterizes the anisotropy of the image, this index varying monotonically as a function of the statistical dispersion of the values of the function τ(θ) for θ varying between 0 and π.

The invention relates to a method for characterizing the anisotropy ofthe texture of a digital image. The invention also relates to a methodfor classifying digital images depending on the anisotropy of theirtexture. The invention lastly relates to a data storage medium and to acomputer for implementing these methods.

Patent application WO2016/042269A1 describes a method that allows theHurst exponent H of the texture of an image and terms β_(j) that varydepending on the characteristics of the texture of this image in oneparticular direction corresponding to an angle α_(j) to be estimated.This method is very good at identifying the anisotropy of an image.

It has also been proposed to construct an index A that characterizes theanisotropy of the texture of the image from the terms β_(i). This indexis called the “anisotropy index”. For example, the calculation of ananisotropy index A from the terms β_(i) is described in the followingarticles:

-   F. J. P Richard: “Analysis of anisotropic brownian textures and    application to lesion detection in mammograms”, Procedia    Environmental Sciences 27 (2015) 16-20, 2015,-   F. J. P Richard: “Anisotropy indices for the characterization of    brownian textures and their application for breast images”, 2016-   F. J. P Richard: “Some anisotropy indices for the characterization    of Brownian texture and their application to breast images”, SPATIAL    STATISTICS, vol. 18, Feb. 17, 2016, pages 147-162,    -   F. J. P Richard: “Tests of isotropy for rough textures of        trended images”, STATISTICA SINICA, vol. 26, no 3, Jul. 1, 2016,        pages 1-32.

In these articles, the index A is dependent on the average of the termsβ_(j). For a given angle α_(j), the term β_(j) varies as a function ofcharacteristics of the texture in this given direction but also as afunction of the Hurst exponent H. It will be recalled that the Hurstexponent H is a global characteristic of the texture that is independentof the orientation of the image. Thus, when a variation in the termβ_(j) is observed, it is not possible to know simply whether thisvariation is due to a modification of the anisotropy of the image or toa modification of the texture in its entirety and therefore of the Hurstexponent H. Thus, this index A varies not only as a function of theanisotropy of the texture of the image but also as a function of theHurst exponent H of the texture of this image.

The invention aims to provide a method for characterizing the anisotropyof an image using an anisotropy index that varies as a function of theanisotropy of the texture while being much less sensitive to variationsin the Hurst exponent H of the same texture. One subject thereof istherefore a method as claimed in claim 1.

The claimed method estimates from the terms β_(j), the coefficients of afunction τ(θ), here referred to as the asymptotic topothesia function.This function τ(θ) has the particularity of returning a value thatcharacterizes the texture of the image in the direction θ while beingalmost completely independent of the value of the Hurst exponent Hassociated with the same texture. Thus, constructing the anisotropyindex, which varies monotonically as a function of the statisticaldispersion of the function τ(θ), in this way allows an index that variesas a function of the anisotropy of the texture while being practicallyindependent of the value of the Hurst exponent H of the same texture, tobe obtained.

Embodiments of this method may have one or more of the features of thedependent claims.

These embodiments of the characterizing method furthermore have thefollowing advantages:

-   -   The calculation of the estimates of the scalar coefficients of        the function τ(θ) by virtue of a linear relationship between        these coefficients to be estimated and the terms β_(j)        accelerates very substantially the execution of the        characterizing method.    -   The fact that the optimal value of the parameter λ is calculated        allows the estimation of the scalar coefficients of the function        τ(θ) to be improved and therefore an index that is even less        sensitive with respect to variations in the Hurst index H to be        obtained.

Another subject of the invention is a method for automaticallyclassifying digital images depending on the anisotropy of their texture.

The invention also relates to a data storage medium, containinginstructions for implementing the claimed method, when theseinstructions are executed by a computer.

The invention also relates to a computer for implementing the claimedmethod.

The invention will be better understood on reading the followingdescription, which is given merely by way of nonlimiting example withreference to the appended drawings, in which:

FIGS. 1A to 1D are schematic illustrations of digital images havingisotropic and anisotropic textures;

FIG. 2 is a schematic illustration of a computing device forautomatically characterizing the anisotropy of a digital image;

FIG. 3 is a flowchart of a method for characterizing the anisotropy ofthe texture of a digital image;

FIG. 4 is a flowchart of a method for automatically classifying imagesdepending on the anisotropy of their texture,

FIGS. 5A to 5F are illustrations of digital images of the texture ofvarious types of paper,

FIG. 6 is a graph showing the distribution of the texture of variouspapers as a function of their Hurst exponent and of their anisotropyindex.

In these figures, the same references have been used to referenceelements that are the same. In the rest of this description, featuresand functions well known to those skilled in the art have not beendescribed in detail.

In this description, the following mathematical notations andconventions have been adopted, unless otherwise indicated:

-   -   the interval [X, Y] is the interval of all the integer numbers        higher than or equal to X and lower than or equal to Y, where X        and Y are themselves integers;    -   a vector A in a space of d dimensions (such that        ^(d)) has for coordinates A₁, A₂, . . . , A_(d);    -   [0, X]^(d) is the product [0,X₁]×[0,X₂]× . . . ×[0,X_(d)] where        X is a vector of        ^(d) of coordinates X₁, X₂, . . . , X_(d), such that the i-th        coordinate U, of a vector U of [0, X]^(d) belongs to the        interval [0,X_(i)], where i is an index higher than or equal to        0 and lower than or equal to d;    -   |X| is the sum of the components of the vector X, such that        |X|=|X₁|+|X₂|+ . . . +|X_(d)|;    -   |X|² is the Euclidean norm to the power of two of the vector X,        such that |X|²=(X₁ ²+X₂ ²+ . . . +X_(d) ²).

FIG. 1A shows a digital image 2 the texture of which has anisotropy.

In this description, anisotropy is understood to mean the fact that theproperties of the texture of the image differ depending on the directionin which they are observed.

The texture of a digital image is generally defined relative to thespatial distribution of variations in the intensity and/or variations inthe tone of the pixels forming the digital image. Texture is amanifestation of the Hölderian regularity of the image. The notion oftexture is for example defined:

-   -   in the work “Handbook of Texture Analysis”, M. Mirmehdi et al.,        eds., World Scientific, October 2008 in the chapter        “Introduction to Texture Analysis” by E. R. Davies, or even:    -   in section “I. Introduction” of the article by Robert M.        Haralick et al; “Textural features for image classification”;        IEEE Transactions on Systems, Man and Cybernetics; vol. SMC-3,        no 6, p. 610-621, November 1973.

The anisotropy of an image may be caused by two factors: texture and“tendency”. Typically, texture corresponds to short-range (i.e.high-frequency) variations in the intensity of the pixels whereastendency relates to longer-range (i.e. low-frequency) variations in theintensity of the pixels.

Here, it is the texture, and above all its anisotropy, that is ofinterest for characterization of the image 2. For example, when theimage 2 represents a biological tissue, the anisotropic character of thetexture of the image may give an indication as to the presence or therisk of development of cancerous cells within this tissue. In thisexample, the image 2 is a mammogram.

FIGS. 1B to 1D illustrate other examples of images that could correspondto a mammogram. FIG. 1B shows an image the texture of which isisotropic. FIGS. 1C and 1D respectively show images the texture of whichis isotropic and anisotropic and that each comprise an anisotropy causedby a polynomial tendency of the second order. This tendency is orientedin the horizontal direction of these images.

Image 2 is formed from a plurality of pixels. Each pixel is associatedwith:

-   -   a pixel intensity value, and with    -   a position p in the space        ^(d)        where d is a natural integer higher than or equal to 2 that        represents the dimension of the image 2. Here, in this example,        d=2.

Thus, spatially the pixels of the image 2 are arranged in a lattice inthe space

^(d). Preferably, the resolution of the image is the same along all thed axes of the image. Below, the set of possible positions of the pixelsof the image 2 is denoted [0, N]^(d), where N is a vector that codes thesize of the image and the components of which are strictly positivenatural integers belonging to

^(d). This notation means that the p₁, p₂, . . . , p_(d) coordinates ofthe position p of a pixel of the image belong, respectively, to the sets[0,N₁], [0, N₂], . . . , [0, N_(d)], where N₁, N₂, . . . , N_(d) are thecoordinates of N. Here, the image 2 has a square shape of (N₁+1)*(N₂+1)size, where N₁+1=N₂+1 and N₁+1 is the length of one side of this squareexpressed in number of pixels. For example, the image 2 is a zone ofinterest extracted from an image of larger size. The sides of the image2 have a length larger than or equal to 50 pixels or to 100 pixels or to500 pixels.

In this example, the light intensity of the pixels is encoded ingrayscale, in 8-bit grayscale for example. The pixel intensity valuesare integer numbers belonging to the interval [0,255].

FIG. 2 shows a device 12 for identifying and characterizing theanisotropy of the texture of the image 2. The device 12 is able, for agiven image 2, to indicate whether the image is isotropic or anisotropicand, advantageously, in the latter case, to quantify, i.e. characterize,the amplitude of the anisotropy.

The device 12 comprises to this end:

-   -   a programmable computer 14, such as a microprocessor,    -   a data storage medium 16, such as a memory,    -   an interface 18 for acquiring a digital image.

The interface 18 allows the image 2 to be acquired. Typically, thedigital image is generated by an electronic image-capturing apparatussuch as a radiography apparatus. The computer 14 executes theinstructions stored in the medium 16. The medium 16 in particularcontains instructions for implementing the method of FIGS. 3 and 4,which will be described in more detail below.

The anisotropy of the image 2 is identified and characterized using acertain number of operations.

In this example, the image 2 is modelled as an intrinsic random gaussianfield. In other words, the intensity value associated with each pixel ofthe image 2 is said to correspond to a Gaussian random variable Z. Theconcept of intrinsic random gaussian fields is defined in more detail inthe following works: J. P. Chilès et al. “Geostatistics: ModelingSpatial Uncertainty”, J. Wiley, 2^(nd) edition, 2012.

The intensity value associated with the pixel the position of which inthe image 2 is given by the position p is denoted Z[p]. For example, anorthonormal coordinate system having as origin the null vector (0)^(d)is defined in

^(d). The position p belongs to

^(d).

An example of an implementation of the method for automaticallycharacterizing the anisotropy of a texture will now be described withreference to the flowchart of FIG. 3 and to FIGS. 1 and 2.

In a step 20, the image 2 is automatically acquired by the interface 18and stored, for example, in the medium 16. The notation “I” will be usedbelow to refer to this image 2.

In this example, of dimension d=2, the normalized image 2 is modelled bya square matrix Z of (N₁+1)*(N₂+1) size. The coefficients of this matrixZ are the Z[p] corresponding to the intensity of the pixels of positionp. The components of the vector p give the position of this coefficientin the matrix Z. For example, Z[p] is the coefficient of the p₁-th rowand the p₂-th column of Z, where p₁ and p₂ are the coordinates of theposition p in [0, N]².

Next, in a step 22, geometric conversions are applied to the image 2 inorder to obtain a series of converted images I_(j,k). These conversionscomprise modifications T_(j,k) of the image 2 that each include:

-   -   a rotation by an angle α_(j), and    -   scaling by a scale factor γ_(k).

Below, the image obtained after the modification T_(j,k) has beenapplied to the acquired image I is denoted T_(j,k)(I).

Each modification T_(j,k) is characterized uniquely by a vector u_(jk)of the space

²\{(0,0)}, such that α_(j)=arg(u_(jk)) and γ_(k)=|u_(jk)|². The space

²\{(0,0)} is the space

² deprived of the point of coordinates (0,0).

The indices “j” and “k” are integer indices that uniquely identify theangle α_(j) and the factor γ_(k), respectively. The index j variesbetween 1 and n_(j). To simplify notation, “rotation j” and “scaling k”will be spoken of below when referring to a rotation of angle α_(j) andto scaling by a factor γ_(k), respectively.

The rotation j makes each of the pixels of the image 2 rotate by theangle α_(j) from a starting position to an end position about the samepredetermined point or the same predetermined axis. Typically, thispoint or this axis of rotation passes through the geometric center ofthe image. The rotation here occurs about the geometric center of theimage. The geometric center of a digital image is defined as being thecentroid of the positions of all of the pixels of the image, eachweighted by a coefficient of the same value.

The scaling k enlarges or shrinks the image via a homothety of factorγ_(k). In the following examples, the center of the homothety is thegeometric center of the image.

These modifications T_(j,k) are applied for at least two, and,preferably, at least three or four angles α_(j) of different values.Advantageously, the various values of the angles α_(j) are distributedas uniformly as possible between 0° and 180° while respecting theconstraint that the vector u_(jk) must belong to the space

²\{(0,0)}. The number n_(j) of different values for the angle α_(j) isgenerally chosen not to be too high in order to limit the number ofcalculations to be carried out. For example, this number n_(j) is chosento be lower than 150 or 100. A good compromise consists in choosing atleast four different values for the angle α_(j) and, preferably at least10 or 20 different values. For each angle α_(j), modifications T_(j,k)are applied for at least two, and, preferably, at least three or four orfive, different scalings γ_(k).

The values of the factor γ_(k) are for example higher than or equal to 1and lower than or equal to 10² or 8² or 4². Preferably, the variousvalues of the factor γ_(k) are distributed as uniformly as possible inthe chosen interval of values. For example, here, the scalings γ_(k)used are all scalings for which the following condition is met: theEuclidean norm of the vector u_(jk) belongs to the interval [√2; 10].

For example, the angles of rotation α_(j) are chosen depending on thehorizontal and vertical directions of the image 2. For example, to maketwo rotations j1 and j2, values α_(j1)=90° and α_(j2)=180° are chosen,where j1 and j2 are particular values of the index j. The angles arehere expressed with respect to the horizontal axis of the image 2.

In this example, of dimension d=2, the modifications T_(j,k) applied arethe following, expressed here in matrix form:

T j , k = k  [ cos  ( α j ) - sin  ( α j ) sin  ( α j ) cos  ( α j) ]

In step 22, K-increments are calculated for each of the converted imagesT_(j,k)(I). This calculation includes filtering intended to removetendencies of polynomial form of order strictly lower than K. Moreprecisely, for each image T_(j,k)(I), a filter is applied allowing theK-increment V_(j,k) of this image T_(j,k)(I) to be calculated. It is theK-increment of this image T_(j,k)(I) that forms the converted imageI_(j,k). The K-increment V_(j,k) of this image is not calculated for allthe points of the image T_(j,k)(I), but only for certain thereof as willbe seen below.

The K-increment concept is for example defined in more detail in thefollowing work: J. P. Chilès et al. “Geostatistics: Modeling SpatialUncertainty”, J. Wiley, 2^(nd) edition, 2012.

Here, the filtering is carried out by means of a convolution kerneldenoted “v”, in order to achieve linear filtering. Below, the “filter v”will be spoken of when referring to this convolution kernel.

The filter v is defined for the set [0,L]^(d). This filter v ischaracterized by a characteristic polynomial Q_(v)(z) defined by:

∀ z ∈ d , Q v  ( z ) = ∑ p ∈ [ 0 , L ] d  v  [ p ]  z p

Here, the filter v is a matrix and the quantity v[p] is a particularscalar value for this filter for the position p, where p is a vector of[0,L]^(d). This value v[p] is zero if the vector p does not belong to[0,L]^(d). Equivalently, the filter v is also said to have a boundedsupport in [0,L]^(d). This filter v is different from the null functionthat, for any value of the vector p, has a value v[p] of zero. Thenotation z^(p) here designates the monomial z₁ ^(p1)*z₂ ^(p2)* . . .*z_(d) ^(pd).

The filter v is thus parameterized with the vector L, which is a vectorof [0, N]^(d). Generally, the vector L is chosen so as to be containedin the image I. Therefore, preferably, values of L are used thatrespect, for any i ranging from 1 to d, the relationship L_(i)<<N_(i),i.e. L_(i) is lower than 10 times or 100 times N_(i).

In addition, the filter v is such that its characteristic polynomialQ_(v)(z) meets the following condition:

∀a∈[0,K]^(d) such that |a|≤K then

${\frac{\partial^{a}Q_{v}}{{\partial z_{1}^{a_{1}}}\mspace{14mu} \ldots \mspace{14mu} {\partial z_{d}^{a_{d}}}}\left( {1,\ldots \mspace{14mu},1} \right)} = 0$

where the constant K is a nonzero natural integer and ∂^(|a|)Q_(v)/∂z₁^(a1) . . . ∂z_(d) ^(ad) is the partial derivative of the polynomialQ_(v)(z) with respect to the components of the vector z, the symbol∂z_(i) ^(ai) indicating a differentiation of the polynomial Q_(v)(z) oforder a_(i) with respect to the variable z_(i), where z_(i) is the i-thcomponent of the vector z and a_(i) is the i-th component of the vectora, i being an integer index higher than or equal to 0 and lower than orequal to d.

Filtering the image T_(j,k)(I) with the filter v allows the effect of“tendency” to be removed from the subsequent calculations of the methodwhen said tendency takes the form of a polynomial of order P_(o),provided that the value of the constant K is chosen as follows:

-   -   K≥P_(o)+1 if d is lower than or equal to 4, and    -   K≥P_(o)/2+d/4 if d>4.

The K-increments of the image T_(j,k)(I), denoted V_(j,k), are thencalculated by virtue of the filter v as follows:

${V_{j,k}\lbrack m\rbrack} = {\sum\limits_{p \in {\lbrack{0,L}\rbrack}^{d}}{{v\lbrack p\rbrack} \cdot {Z\left\lbrack {m - {T_{j,k} \cdot p}} \right\rbrack}}}$

where:

-   -   V_(j,k)[m] is a K-increment calculated in the image T_(j,k)(I)        for the pixel of position m, with m a vector belonging to a set        E that will be defined below;    -   the product T_(j,k)·p corresponds to the application of the        modification T_(j,k) to the pixel of position p of the image I        and expresses the coordinates in        ^(d), after application of the modification T_(j,k), of the        pixel that initially had the position p in the image I,    -   v[p] is the value of the filter v for the value of p.

For each image T_(j,k)(I), the K-increment is calculated only for thosepixels of the image T_(j,k)(I) the positions of which belong to a set E.The set E only contains positions:

-   -   that belong to the image I, and    -   that, whatever the applied modification T_(j,k), occupy a        position that already exists in the image I, after application        of this modification T_(j,k).

Furthermore, for any position m belonging to E, the pixels of positions“m−T_(j,k)·p” occupy a position contained in the image I.

The number of positions that belong to the set E is denoted n_(E).

Thus, quadratic variations are calculated solely for those points of theconverted image for which no interpolation is necessary. It is thuspossible to make use of rotations j of any angle, in contrast to thecase of projections. Specifically, if a projection is carried out in a,for example, diagonal direction of the image, some projected points havea position that does not belong to the set [0,N]^(d). In other words,they do not form part of the lattice. It is therefore necessary toemploy an interpolation in order to determine the intensity valueassociated with the points that do belong to this lattice. Thisintroduces an approximation and therefore an error. With themodifications T_(j,k) and then the selection of the points of the set E,the reliability of the method is improved.

In this example, the filtering is carried out within the same formula asthe application of the modifications T_(j,k).

With this choice of the constant K, the filtering generates incrementsV_(j,k)[m] of order K. This filtering makes it possible to not take intoaccount any anisotropy of the images caused by tendency, but rather onlythe anisotropy of the texture of the underlying image. As a result, thereliability of the characterizing method is improved.

The constant K must therefore be chosen as described above depending onthe nature of the tendency present in the image 2. Typically, in thecase of a mammogram, the polynomial degree P_(o) of the tendency islower than or equal to 2. For example, a value is selected by a user ofthe method. To this end, step 22 here comprises acquiring a value of thevector L and a value of the constant K.

In this example, the filter v is chosen as follows:

${v\lbrack p\rbrack} = {{\left( {- 1} \right)^{p} \times C_{L_{1\;}}^{p_{1}} \times \ldots \times C_{L_{d}}^{p_{d}}} = {\left( {- 1} \right)^{p} \times {\prod\limits_{i = 1}^{d}\frac{L_{i}!}{{p_{i}!} \times {\left( {L_{i} - p_{i}} \right)!}}}}}$

if the vector p belongs to [0,L]^(d) and v[p]=0 if not, where the termsC^(p) _(L) are binomial coefficients.

With this particular filter, the condition expressed above on thecharacteristic polynomial Q_(v)(z) is met if K=|L|−1. Thus, the value ofK is deduced from the acquired value of the parameter L.

Thus, in this particular case, filtering the image T_(j,k)(I) with thefilter v allows the effect of “tendency” to be removed when the latterhas a polynomial form of order P_(o), provided that the parameter L ischosen as follows:

-   -   |L|=P_(o)+2 if d is lower than or equal to 4, and    -   |L|=P_(o)/2+d/4+1 if d>4.

In this example, of dimension d=2, the vector L has two components, L₁and L₂. To remove a tendency of polynomial degree P_(o)=2, it isnecessary to choose L₁ and L₂ such that |L| is equal to four. PreferablyL₁ and L₂ are chosen such that L₁=4 and L₂=0. Specifically, by choosingvalues, for the coordinates of the vector “L”, that are sufficientlydifferent from each other, the directional sensitivity of the filter isincreased. Thus, it will react more markedly, and will therefore filtermore effectively, variations that are oriented in one particulardirection. In contrast, a filter for which L₁ and L₂ are chosen suchthat L₁=2 and L₂=2 would be less sensitive to directional signals andless effective.

In this embodiment, the filter v is defined by the followingrelationship:

$\begin{matrix}{{v\lbrack p\rbrack} = {\left( {- 1} \right)^{p_{i}}\frac{L_{1}!}{{p_{1}!}{\left( {L_{1} - p_{1}} \right)!}}}} & (1)\end{matrix}$

if the vector p belongs to [0,L]×{0} and v[p]=0 if not, where L₁ belongsto

\{0}. Here, L₁=4. Under these conditions, the order of the kernel v isequal to K=L₁−1.

In this step, step 22, for each different value of j and k, the computer14 carries out in succession the following operations:

-   -   application of the modification T_(j,k) to the image 2 in order        to obtain the image T_(j,k)(I), then    -   application of the filter v to the image T_(j,k)(I) in order to        obtain the converted image I_(j,k).

Next, in a step 24, for each image I_(j,k), the p-variation W_(j,k)associated with this image I_(j,k) is calculated. The concept ofp-variations is well known to those skilled in the art of statistics andprobabilities. Here, it is calculated in the following way:

${W_{j,k}\lbrack m\rbrack} = {\frac{1}{n_{E\;}}{\sum\limits_{m \in E}\left( {V_{j,k}\lbrack m\rbrack} \right)^{q}}}$

In the above equation, the symbol “q” is used instead of the symbol “p”conventionally used in this equation in order to avoid any confusionwith the symbol “p” used in this description to designate the positionof a pixel. In this example, a particular form of p-variations is used:“quadratic variations” or “2-variations”, for which q=2. Thus, thequadratic variation W_(j,k) of the image I_(j,k) is calculated in thefollowing way, from the K-increments calculated after filtering in step22:

${W_{j,k}\lbrack m\rbrack} = {\frac{1}{n_{E\;}}{\sum\limits_{m \in E}\left( {V_{j,k}\lbrack m\rbrack} \right)^{2}}}$

These quadratic variations W_(j,k) contain information that is importantfor the identification of anisotropy. To extract this information, theprocedure is as follows.

In a step 26, a covariance analysis comprising a statistical regressionis carried out on all the variations W_(j,k) calculated for each of theimages I_(j,k) in order to estimate:

-   -   the value of the Hurst exponent H of the image I and,    -   a term β_(j).

The statistical regression is defined by the following relationship:

log(|W _(j,k)|)=log(|u _(jk)|²)*H+β _(j)+ε_(j,k),

where:

-   -   |u_(jk)|² is the Euclidean norm to the power of two of the        vector u_(jk),    -   H is the Hurst exponent of the image I;    -   β_(j) is a quantity that does not depend on the scaling k; here        this parameter is analogous to a so-called intercept parameter        of the regression, unless it depends on the rotations j;    -   ε_(j,k) is an error term of the regression the statistical        properties of which are predetermined and set by the user. For        example, the error terms ε_(j,k) are intercorrelated Gaussian        random variables.

It will be recalled that the Hurst exponent H is a physical quantitythat is independent of the rotations of the image.

Thus, a number n_(j) of terms β_(j) is obtained, n_(j) being the numberof different rotations applied to the image I.

For example, if it is decided to make do with only the two rotations j1and j2 described above, the regression is carried out on the basis ofall the quadratic variations calculated for j1 and for j2. Thus, twoterms β_(j1) and β_(j2) are obtained.

At this stage, in a step 28, the computer 14 estimates the scalarcoefficients τ_(m) of an even function τ(θ) called the asymptotictopothesia function. The function τ(θ) is continuous in the interval [0;2π]. The function τ(θ) is a function that minimizes the followingcriterion C:

$\begin{matrix}{C = {\sum\limits_{j = 1}^{n_{i}}\left( {\beta_{j} - {\tau*{\Gamma \left( \alpha_{j} \right)}}} \right)^{2}}} & (2)\end{matrix}$

where:

-   -   the symbol “τ*Γ(θ)” designates the function obtained by carrying        out the circular convolution of the functions τ(θ) and Γ(θ),    -   Γ(θ) is the function defined by the following relationship:

Γ(θ) = ∫_(ℝ+)v̂(ρ θ)²ρ^(−2H − 1)d ρ

where:

-   -   {circumflex over (V)} is the discrete Fourier transform of the        kernel v,    -   H is the Hurst exponent of the acquired image,    -   ρ is the variable of integration.

The scalar coefficients τ_(m) of the function τ(θ) are defined by thefollowing relationship:

${\tau (\theta)} = {\sum\limits_{m = 0}^{M}{\tau_{m}{f_{m}(\theta)}}}$

where:

-   -   M is an integer number higher than one,    -   τ_(m) are the scalar coefficients of the function τ(θ),    -   f_(m)(θ) are the functions of a basis of π-periodic functions        defined in the interval [0; 2π].

A π-periodic function is a periodic function of period π.

In this embodiment, the π-periodic function basis used is a Fourierbasis. Therefore, the function τ(θ) is here defined by the followingrelationship:

${\tau (\theta)} = {\tau_{0} + {\sum\limits_{m = 1}^{M}\left( {{\tau_{1,m}{\cos \left( {2m\; \theta} \right)}} + {\tau_{2,m}{\sin \left( {2\; m\; \theta} \right)}}} \right)}}$

where τ₀, τ_(1,m) and τ_(2,m) are the scalar coefficients of thefunction τ(θ),

The number M is a constant that is defined beforehand, for example bythe user. Generally, this number M is lower than the number n_(j) ofdifferent angles α_(j). Typically, this number M is also higher than orequal to 2 or 4. Typically, the number M is chosen so that the number ofscalar coefficients of the function τ(θ) is comprised in the interval[0.35n_(j), 0.75n_(j)] or in the interval [0.45n_(j), 0.55n_(j)].

In this embodiment, and in the context specified above, it has beenpossible to establish a linear relationship between an approximation τ*of the coefficients of the function τ(θ) and the estimated terms β_(j).This relationship is the following:

τ*=(L ^(T) L+λR)⁻¹ L ^(T) β  (3)

where:

-   -   τ* is the vector (τ₀*, τ_(1,1)*, τ_(2,1)*, τ_(1,2)*, τ_(2,2)*, .        . . , τ_(1,M-1)*, τ_(2,M-1)*, τ_(1,M)*, τ_(2,M)*)^(T), the        coefficients τ₀*, τ_(1,1)*, τ_(2,1)*, τ_(1,2)*, τ_(2,2)*, . . .        , τ_(1,M-1)*, τ_(2,M-1)*, τ_(1,M)*, τ_(2,M)* being the estimates        of the coefficients τ₀, τ_(1,1), τ_(2,1), τ_(1,2), τ_(2,2), . .        . , τ_(1,M-1), τ_(2,M-1), τ_(1,M), τ_(2,M), respectively    -   L is a matrix of n_(j)×(2M+1) size the k-th column of which is        defined by the following relationship:

({circumflex over (μ)}[0],{circumflex over (μ)}[1]cos(α_(k)),{circumflexover (μ)}[1]sin(α_(k)), . . . ,{circumflex over (μ)}[M]cos(Mα_(k)),{circumflex over (μ)}[M] sin(Mα _(k)))

-   -   -   where {circumflex over (μ)}[0], {circumflex over (μ)}[1], .            . . , μ[M] are the coefficients of the discrete Fourier            transform of the following function μ(θ):            μ(θ)=|cos(θ)|^(2H), where H is the Hurst exponent of the            acquired image,

    -   λ is a predetermined parameter,

    -   R is a diagonal matrix of (2M+1)×(2M+1) size the coefficients of        which on the diagonal are, in order: 0, 2, 2, 5, 5, . . . ,        (1+M²), (1+M²),

    -   the symbol “^(T)” designates the transpose operation,

    -   β is the vector (β₁, β₂, . . . , β_(nj-1), β_(nj))^(T).

The inventors have also established that the optimum value of theparameter λ is equal or very close to a value λ*. The value λ* isobtained using the following relationship:

$\begin{matrix}{\lambda^{*} = {\frac{\kappa}{\left( {1 + M^{2}} \right)}\frac{{trace}\left( {V\left( \underset{\_}{\beta} \right)} \right)}{{\underset{\_}{\beta}}^{2}}}} & (4)\end{matrix}$

where:

-   -   κ is equal to v₊/v⁻, v₊ and v⁻ being the lowest and the highest        eigenvalues of the matrix L^(T)L, respectively,    -   trace(X) is a function that returns the sum of the diagonal        coefficients of a square matrix X,    -   V(β) is the covariance matrix of the estimator of the vector β,    -   |β|² is the Euclidean norm to the power of two of the vector β.

Thus, in this embodiment, in step 28, the computer 14 calculates thevalue λ* using the above relationship. Next, it chooses a value of theparameter λ close to the value λ*. For example, it chooses, randomly forexample, a value of the parameter λ in the interval [0; 1.3λ*] or [0;1.1λ*]. Most often, the value of the parameter λ is chosen in theinterval [0.7λ*; 1.3λ*] or [0.9λ*; 1.1λ*]. Here, the value of theparameter λ is systematically set equal to the value λ*. Lastly, thecomputer estimates the coefficients τ₀, τ_(1,m), τ_(2,m) usingrelationship (3).

The relationship (4) was established by searching for a numericalexpression for the coefficients τ₀*, τ_(1,m)*, τ_(2,m)* that minimizesnot the criterion C directly but a penalized criterion C_(A). Thispenalized criterion C_(A) is for example the following:

$\begin{matrix}{C_{A} = {C - {\lambda {\sum\limits_{m = 1}^{M}\; {\left( {1 + m^{2}} \right)\left( {\tau_{1,m}^{2} + \tau_{2,m}^{2}} \right)}}}}} & (5)\end{matrix}$

The term that is subtracted from the criterion C in relationship (4) isknown as the “penalty”.

In addition, in the case where the filter v is that defined byrelationship (1), then the circular convolution τ*Γ(θ) may be written inthe following form:

τ*Γ(θ)=

τ*μ(θ)  (6)

where:

-   -   μ(θ) is defined by the following relationship:

μ(θ)=|cos(θ)|^(2H)

-   -   γ is a constant that is independent of the value θ,    -   the symbol “*” is a circular convolution such that τ*μ(θ) is the        function resulting from the circular convolution of the        functions τ(θ) and μ(θ).

The constant γ is defined by the following relationship:

 = 2∫_(ℝ+)v̂(ρ)²ρ^(−2H − 1)dρ

where the symbols used in the above relationship were defined above and| . . . |² is the Euclidean norm to the power of two.

Thus, the penalized criterion C_(A) may be written in the followingform:

C _(A) =|Lτ−β| ²+λτ^(T) Rτ

where the symbols L, τ, β, λ and R were defined above and | . . . |² isthe Euclidean norm to the power of two. The vector τ* is the numericalsolution of the criterion C_(A).

For a given value of the angle θ comprised between zero and π, the valueof the function τ(θ) depends on the characteristics of the texture ofthe image in the direction θ. This value is independent of the value ofthe Hurst exponent H. Thus, the statistical dispersion of the values ofthe function τ(θ) for θ varying between 0 and π is representative of theanisotropy of the texture of the image. For example, the statisticaldispersion of the function τ(θ) is represented by an index A that isdependent on the sum of the following deviations: |τ−M_(τ)|_(Lp) for θvarying between 0 and π, where:

-   -   τ is the function τ(θ),    -   M_(τ) is an average value of the values of the function τ(θ) for        θ varying from 0 to π or an approximation of this average, and    -   | . . . |_(Lp) is the norm L1 if Lp is equal to 1, the norm L2        if Lp is equal to 2 and so on, Lp being strictly higher than        zero.

Thus, the statistical dispersion of the function τ(θ) may be thevariance or the standard deviation of the values of this function in [0;π]. For example, here, Lp is chosen equal to 2 and the calculated indexA is equal to the square root of the sum defined above. Thus, the higherthe value of the index A, the greater the anisotropy of the image. Underthese conditions, the index A is defined by the following relationship:

$\begin{matrix}{A = \sqrt{\int_{\lbrack{0,{\pi\lbrack}}}{\left( {{\tau (ɛ)} - M_{\tau}} \right)^{2}{dɛ}}}} & (7)\end{matrix}$

Here, in a step 30, the computer 14 calculates the index A using thefollowing formula, which corresponds to relationship (7):

$A = \sqrt{\sum\limits_{m = 1}^{M}\; \left( {\tau_{1,m}^{\;^{*}2} + \tau_{2,m}^{\;^{*}2}} \right)}$

The flowchart of FIG. 4 describes an example of application of the abovemethod to the automatic classification of images 2 with respect to oneanother depending on their texture. This example is given in theparticular case where the images are photographs of sheets of papertaken using a microscope and under grazing illumination. Such capturedimages are illustrated in FIGS. 5A to 5F. FIGS. 5A to 5C show images ofsheets of gloss paper from three different manufacturers. FIGS. 5D and5E show images of sheets of satin paper. FIG. 5F shows an image of asheet of matte paper. The databases containing these images aredescribed in detail in the following articles:

-   R. Johnson, P. Messier, W. Sethares, et al: “Pursuing automated    classification of historic photographic papers from raking light    images”, J. AM. Inst. Conserv. 53(3):159-170, 2014, and-   P. Messier, R. Johnson, H. Wilhelm, W. Sethares, A. Klein, et al:    “Automated surface texture classification of inkjet and photographic    media”, In NIP & Digital Fabrication Conference, pages 85-91,    Society for Imaging Science and Technology, 2013.

In a step 40, a plurality of digital images 2 are automaticallyacquired. Among the acquired images, certain correspond to gloss paper,others to satin paper and others to matte paper.

In a step 42, for each of these images, the anisotropy index A and theHurst exponent H are calculated by implementing the method of FIG. 3. Inthis particular case, the method of FIG. 3 was implemented with M=23 andn_(j)=96.

In a step 44, the acquired images are automatically classified withrespect to one another depending on their index A and exponent Hcalculated in step 42. This classification is for example achieved bymeans of a classifier such as a classifying algorithm based on neuralnetworks or a support vector machine.

The graph of FIG. 6 shows for each image a point of coordinates (H, A),where H and A are the Hurst exponent and anisotropy index calculated forthis image, respectively. Here, the function τ(θ) has been normalized.In this graph, the points represented by crosses, circles and rhombicorrespond to images of a piece of gloss, satin and matte paper,respectively. This graph shows that the combination of the Hurstexponent H and of the anisotropy index A allows the various types ofpaper to be effectively distinguished between. Here, all the pieces ofgloss, satin and matte paper are in very different zones. These zonesare encircled in the graph of FIG. 6. In addition, inside a given zone,a bunch of points grouped together in an even smaller zone oftencorresponds to one particular manufacturer or to one particular batch.Thus, the classification is able not only to distinguish between varioustypes of paper but also, for a given type of paper, distinguish betweendifferent manufacturers or different batches.

Variants of the Image:

The pixels of the image 2 may have other intensity values. The intensityvalue of each pixel may be a real value. Or it may be higher than 256.For example, the image 2 is color-coded. In this case, the color imageis separated into a plurality of monochromatic images each correspondingto one of the color channels used to form the color image. The method isthen applied separately to each of these monochromatic images.

The image 2 may not have a square shape. For example, when of dimensiond=2, the image may have a rectangular or even trapezoidal shape. Whenthe image does not have a regular shape, the notions of “horizontal” and“vertical” direction are replaced by reference directions suited to thegeometry of the image. For example, in the case of an image oftriangular shape, the base and height of the triangle will be used asreference.

The dimension d of the images may be higher than two. For example, theimage 2 may be a hypercube of dimension d.

The image 2 may be other than a mammogram or a photograph or a sheet ofpaper. For example, it may be a question of a photograph of a bonetissue. The anisotropy of the texture of the image then providesinformation on the presence of bone pathologies, such as osteoporosis.Other broader fields of application may be envisioned, such as othertypes of biological tissues, aerial or satellite images, geologicalimages, or photographs of materials. Generally, the method is applicableto any type of irregular and textured image such as an image obtainedfrom any image-capturing electronic apparatus.

Variants of the Method for Estimating the Terms β_(j) and H:

Other modifications T_(j,k) may be used. For example, when the dimensiond=3, the modifications T_(j,k) cause a rotation j about a given axis ofrotation and a scaling k in a given direction of the image. For example,the following modifications may be used:

T j , k = k  [ cos  ( α j ) - sin  ( α j ) 0 sin  ( α j ) cos  ( αj ) 0 0 0 y k ] T j , k = k  [ cos  ( α j ) 0 - sin  ( α j ) 0 y k 0sin  ( α j ) 0 cos  ( α j ) ] T j , k = k  [ y k 0 0 0 cos  ( α j) - sin  ( α j ) 0 sin  ( α j ) cos  ( α j ) ]

-   -   The above modifications T_(j,k) cause a rotation about an axis        and a scaling in a direction that is not parallel to this axis.

The values of the angle α_(j) may be different. Preferably, values ofthe angle α_(j) that do not require interpolations are chosen. However,it is also possible to choose values of the angle α_(j) such that it isnecessary to interpolate the pixels of the converted image to find thevalues associated with each position p comprised in the set E.

As a variant, the rotation and the scaling are not applied at the sametime.

Other filters v may be used to calculate the K-increments. For example,to calculate the K-increments it is also possible to use any filter theconvolution kernel of which is defined as being equal to theconvolution:

-   -   of any convolution kernel v1, and    -   of a convolution kernel v2 equal to the kernel v described        above.        In the particular case where the kernel v1 is an identity        matrix, the filter v described above is obtained. In contrast,        choosing a kernel v1 different from the identity matrix makes it        possible to construct many different filters from the filters v        described above, all suitable for calculating the K-increments.

The filtering may be implemented differently in step 22. In particular,the conversion and the filtering are not necessarily appliedsimultaneously, but in separate formulas.

As a variant, all the conversions T_(j,k) are first applied to the imageI, then, subsequently, the filters are applied to each of the imagesT_(j,k)(I).

The value of K may be different. In particular, with the filter v chosenin the example, when the image I has no tendency, i.e. when M equals 0,then preferably |L|=2 or, if d>4, |L|=1+d/4.

As a variant, a plurality of filters v_(i) are applied to each imageT_(j,k)(I) in step 22. The number of different filters v_(i) applied toa given image T_(j,k)(I) is denoted R. In this case, the converted imageobtained by applying filter v_(i) to image T_(j,k)(I) is denotedI_(i,j,k) and the K-increment of this image at position m in this imageis denoted V_(i,j,k)[m] where “i” is an index that uniquely identifiesthe applied filter v_(i). This index “i” is here different from theindex “i” used above as dummy variable in particular with reference tothe partial derivative of the polynomial Q_(v)(z). Therefore, it is thuspossible to calculate a plurality of quadratic variations for each imageT_(j,k)(I), one for each filter v_(i) applied to this image T_(j,k)(I).Thus, the quadratic variation calculated for the image I_(i,j,k) isdenoted W_(i,j,k).

To this end, step 22 comprises an operation of selecting filters v_(i),for example from a predefined filter library.

The quadratic variations W_(i,j,k) are then calculated in step 24 usingthe following relationship:

${W_{i,j,k}\lbrack m\rbrack} = {\frac{1}{n_{E}}{\sum\limits_{m \in E}\left( {V_{i,j,k}\lbrack m\rbrack} \right)^{2}}}$

In step 26, the regression is then carried out in the following way, onaccount of the n_(i) applied filters:log(|W_(i,j,k)|)=log(|μ_(j,k)|²)*H+β_(i,j)+ε_(i,j,k), where:

-   -   β_(i,j) is the term β_(j) associated with the filter v_(i);    -   ε_(i,j,k) is the error term ε_(j,k) associated with the filter        v_(i).

Thus, a number n_(b) of terms β_(i,j) is obtained, wheren_(b)=n_(j)*n_(i), n_(j) being the number of different rotations appliedto the image I. In step 28, the method of FIG. 3 is implemented for eachfilter v_(i). Thus, scalar coefficients are obtained for i functionsτ_(i)(θ). The anisotropy index A is then calculated from theapproximated coefficients of each of these functions τ_(i)(θ). Forexample, in a simplified embodiment, an anisotropy index A_(i) iscalculated as described above for each of the functions τ_(i)(θ), thecalculated index A then being the average of these indices A_(i).

When a plurality of filters v_(i) are used, the number of filters v_(i)applied may vary from one image T_(j,k)(I) to the next, provided howeverthat to one filter i there correspond at least two rotations j and, foreach of these rotations j, at least two scalings k.

Variants of the Estimation of τ*:

The penalty used in the criterion C_(A) may be different. Provided thatthe penalty is a differentiable function, then it is possible todetermine a linear relationship, such as relationship (3), that directlyexpresses the estimation of the coefficients τ_(m) as a function of theterms β_(j). In particular, it is possible to find such a linearrelationship whatever the filter v and the basis of π-periodic functionsf_(m)(θ) used. The filter v may therefore be different from that definedby relationship (1) and the basis used may also be different from theFourier basis. When the filter v is different from that defined byrelationship (1) or when the basis is different from the Fourier basis,the linear relationship is different from that defined by relationship(3).

The penalty used in the criterion C_(A) may also be a non-differentiablefunction. In this case, it may be difficult or even impossible toestablish a linear relationship between the estimate of the coefficientsτ_(m) and the terms β_(j). For example, the penalty may use a norm L1 ofthe function τ(θ) that is non-differentiable. In this case, othermethods may be used to approximate the coefficients of the function τ(θ)that minimizes this penalized criterion. For example, the estimate ofthe coefficients τ₀, τ_(1,m), τ_(2,m) that minimize the criterion C_(A)are estimated by executing a known minimization algorithm such as theiterative shrinkage-thresholding algorithm (ISTA) or the fast iterativeshrinkage-thresholding algorithm (FISTA).

The variant described here allows the values of the coefficients τ₀,τ_(1,m), τ_(2,m) that minimize the criterion C_(A) to be estimatedwithout requiring to do so a linear numerical relationship, such asrelationship (3), that allows an estimate of these coefficients T₀,T_(1,m), τ_(2,m) to be directly obtained from the values of the termsβ_(j).

The above embodiment was described in the particular case where thefunction τ(θ) is decomposed in a Fourier basis of π-periodic functions.However, the method described above also works if the function τ(θ) isdecomposed in any other basis of π-periodic functions f_(m) in whicheach function f_(m) is defined in the interval [0; 2π]. In particular,it is not necessary for the basis of the π-periodic functions f_(m) tobe an orthogonal basis. Thus, in the basis of the f_(m) functions, thefunction τ(θ) is defined by the following relationship:

${\tau (\theta)} = {\sum\limits_{m = 0}^{M}\; {\tau_{m}{f_{m}(\theta)}}}$

where the T_(m) are the scalar coefficients of the function τ(θ). Forexample, as a variant, the functions f_(m) are piecewise constantπ-periodic functions in [0; π]. A piecewise constant function is afunction that takes constant values in a plurality of immediatelysuccessive sub-intervals comprised between [0; π].

The method for minimizing the criterion C or the criterion C_(A) bymeans of known algorithms for minimizing such a criterion may beimplemented whatever the form of the adopted function f_(m).

As a variant, the number M may be higher than or equal to the numbern_(j).

Variants of the Calculation of the Index A:

As a variant, the index A is calculated only for angles θ equal to theangles α_(j) and not for all the values of θ comprised between 0 and π.In this case, for example, the index A is only dependent on the sum ofthe following deviations:

$\sum\limits_{j = 1}^{n_{j}}\; {{{\tau \left( \alpha_{j} \right)} - M_{\tau}}}^{no}$

The classification may be carried out differently in step 42. Forexample, the order of classification of the images may be chosendifferently.

1. A method comprising characterizing the anisotropy of the texture of adigital image, wherein characterizing comprises: a) acquiring a digitalimage formed from pixels, each pixel being associated with a lightintensity and with a position in space

^(d), where d is a natural integer higher than or equal to two; b)automatically converting) the acquired image in order to obtain aconverted image I_(j,k), the conversion comprising applying amodification T_(j,k) of the image that makes each pixel of the acquiredimage rotate from one position to the next by an angle α_(j) about apoint or an axis and that enlarges or shrinks the image by a factorγ_(k), where α_(j)=arg(u_(jk)) and γ_(k)=|u_(jk)|², u_(jk) being avector that completely characterizes the modification T_(j,k), theindices j and k uniquely identifying the angle α_(j) and the factorγ_(k), respectively, then, for each converted image, calculating aK-increment V_(j,k)[m] for each pixel of position m of the convertedimage, this K-increment being calculated by applying a convolutionkernel v by means of the following formula:${V_{j,k}\lbrack m\rbrack} = {\sum\limits_{p \in {\lbrack{0,L}\rbrack}^{d}}{{v\lbrack p\rbrack} \cdot {Z\left\lbrack {m - {T_{j,k} \cdot p}} \right\rbrack}}}$where: the product T_(j,k)·p corresponds to the application of themodification T_(j,k) to the pixel that initially had the position p inthe image I; the convolution kernel v performs linear filtering andpossesses a characteristic polynomial Q_(v)(z) and a finite support[0,L]^(d), v[p] being the value of the convolution kernel v for theposition p, the characteristic polynomial Q_(v)(z) being defined by thefollowing formula: ∀ z ∈ d , Q v  ( z ) = ∑ p ∈ [ 0 , L ] d   v  [ p]  z p and meeting the following condition: ∀a∈[0,K]^(d) such that|a|≤K then${\frac{\partial^{a}Q_{v}}{{\partial z_{1}^{a_{1}}}\mspace{14mu} \ldots \mspace{14mu} {\partial z_{d}^{a_{d}}}}\left( {1,\ldots \mspace{14mu},1} \right)} = 0$where: L is an acquired vector of [0, N]^(d) that parametrizes thekernel v, N is a vector belonging to

^(d) that codes the size of the image and the components of which arestrictly positive natural integers; the constant K is a nonzero acquirednatural integer; z is a vector of components z₁, z₂, . . . , z_(d);z^(p) designates the monomial z₁ ^(p1)*z₂ ^(p2)* . . . *z_(d) ^(pd);∂^(|a|)Q_(v)/∂z₁ ^(a1) . . . ∂z_(d) ^(ad) is the partial derivative ofthe polynomial Q_(v)(z) with respect to the components of the vector z,the symbol ∂z_(i) ^(ai) indicating a differentiation of the polynomialQ_(v)(z) of order a_(i) with respect to the variable z_(i), where z_(i)is the i-th component of the vector z and a_(i) is the i-th component ofthe vector a, i being an integer index higher than or equal to 0 andlower than or equal to d; the step b) being executed with n_(j)different angles α_(j) and, for each angle α_(j), with at least twodifferent factors γ_(k), n_(j) being an integer higher than or equal to2 so as to obtain at least four different converted images I_(j,k); c)for each different converted image I_(j,k), calculating a p-variationW_(j,k) of this converted image from said calculated K-increments; d)estimating the terms β_(j) of the following statistical regression:log(|W_(j,k)|)=log(|u_(jk)|²)*H+β_(j)+ε_(j,k), where: H is the Hurstexponent of the acquired image; ε_(j,k) is an error term of theregression the statistical properties of which are predetermined;wherein the method further comprises: e) estimating the scalarcoefficients τ_(m) of an even function τ(θ) defined in [0; 2π] thatminimizes the following criterion C:$C = {\sum\limits_{j = 1}^{n_{j}}\; \left( {\beta_{j} - {\tau*{\Gamma \left( \alpha_{j} \right)}}} \right)^{2}}$where: β_(j) are the terms estimated in step d), τ(θ) is the functiondefined by the following relationship for any angle θ belonging to [0;2π]:${\tau (\theta)} = {\sum\limits_{m = 0}^{M}\; {\tau_{m}{f_{m}(\theta)}}}$where: M is a constant integer number higher than an acquired datum,τ_(m) are the scalar coefficients of the function τ(θ), f_(m)(θ) are thefunctions of a basis of π-periodic functions defined in the interval [0;2π], Γ(θ) is the function defined by the following relationship:Γ(θ) = ∫_(ℝ+)v̂(ρθ)²ρ^(−2H − 1)dρ where: {circumflex over (V)}is the discrete Fourier transform of the kernel v, H is the Hurstexponent of the acquired image, ρ is the variable of integration, thesymbol “*” designates the circular convolution of the functions τ(θ) andΓ(θ), f) then, calculating, depending on the estimate of the scalarcoefficients τ_(m), an anisotropy index that characterizes theanisotropy of the image, this index varying monotonically as a functionof the statistical dispersion of the values of the function τ(θ) for θvarying between 0 and π.
 2. The method as claimed in claim 1, whereinthe kernel v used in step b) is equal to the convolution of anyconvolution kernel and the kernel defined in the following way:${v\lbrack p\rbrack} = {{\left( {- 1} \right)^{p} \times C_{L_{1}}^{p_{1}} \times \ldots \times C_{L_{d}}^{p_{d}}} = {\left( {- 1} \right)^{p} \times {\prod\limits_{i = 1}^{d}\frac{L_{i}!}{{p_{i}!} \times {\left( {L_{i} - p_{i}} \right)!}}}}}$if the vector p belongs to [0,L]^(d) and v[p]=0 if not, where the termsC^(p) _(L) are binomial coefficients, the constant K then being equal toK=|L|−1.
 3. The method as claimed in claim 2, wherein: the functionsf_(m) are functions of the Fourier basis and the function τ(θ) isdefined by the following relationship:${\tau (\theta)} = {\tau_{0} + {\sum\limits_{m = 1}^{M}\left( {{\tau_{1,m}{\cos \left( {2m\; \theta} \right)}} + {\tau_{2,m}{\sin \left( {2m\; \theta} \right)}}} \right)}}$where τ₀, τ_(1,m) and τ_(2,m) are the scalar coefficients of thefunction τ(θ), the kernel v is defined by the following relationship:${v\lbrack p\rbrack} = {\left( {- 1} \right)^{p_{1}}\frac{L_{1}!}{{p_{1}!}{\left( {L_{1} - p_{1}} \right)!}}}$in step e), the estimate of these coefficients is calculated using thefollowing relationship:τ*=(L ^(T) L+λR)⁻¹ L ^(T) β where: τ* is the vector (τ₀*, τ_(1,1)*,τ_(2,1)*, τ_(1,2)*, τ_(2,2)*, . . . , τ_(1,M-1)*, τ_(2,M-1)*, τ_(1,M)*,τ_(2,M)*)^(T), the coefficients τ₀*, τ_(1,1)*, τ_(2,1)*, τ_(1,2)*,τ_(2,2)*, . . . , τ_(1,M-1)*, τ_(2,M-1)*, τ_(1,M)*, τ_(2,M)* being theestimates of the coefficients τ₀, τ_(1,1), τ_(2,1), τ_(1,2), τ_(2,2), .. . , τ_(1,M-1), τ_(2,M-1), τ_(1,M), τ_(2,M), respectively, L is thematrix of n_(j)×(2M+1) size the k-th column of which is defined by thefollowing relationship:({circumflex over (μ)}[0],{circumflex over (μ)}[1]cos(α_(k)),{circumflexover (μ)}[1]sin(α_(k)), . . . ,{circumflex over (μ)}[M]cos(Mα_(k)),{circumflex over (μ)}[M] sin(Mα _(k))) where {circumflex over(μ)}[0], {circumflex over (μ)}[1], . . . , μ[M] are the coefficients ofthe discrete Fourier transform of the function μ(θ) given by:μ(θ)=|cos(θ)|^(2H), where H is the Hurst exponent of the acquired image,λ is a predetermined parameter, R is a diagonal matrix of (2M+1)×(2M+1)size the coefficients of which on the diagonal are, in order: 0, 2, 2,5, 5, . . . , (1+M²), (1+M²), the symbol “^(T)” designates the transposeoperation, β is the vector (β₁, β₂, . . . , β_(nj-1), β_(nj))^(T). 4.The method as claimed in claim 3, wherein the method comprises:calculating a value λ* defined by the following relationship:$\lambda^{*} = {\frac{\kappa}{\left( {1 + M^{2}} \right)}\frac{{trace}\; \left( {V\left( \underset{\_}{\beta} \right)} \right)}{{\underset{\_}{\beta}}^{2}}}$wherein κ is equal to v₊/v⁻, v₊ and v⁻ being the highest and the lowesteigenvalue of the matrix L^(T)L, respectively trace(X) is the functionthat returns the sum of the diagonal coefficients of a square matrix X,V(β) is the covariance matrix of the vector β, |β|² is the Euclideannorm to the power of two of the vector β, and automatically choosing theparameter λ in the interval [0; 1.3λ*].
 5. The method as claimed inclaim 1, wherein the anisotropy index is calculated from the sum of thefollowing deviations:$\sum\limits_{i = 1}^{n_{i}}{{{\tau \left( \alpha_{j} \right)} - M_{\tau}}}_{Lp}$wherein M_(τ) is an estimate of an average value of the function τ(θ)for θ varying between 0 and π, | . . . |_(Lp) is the norm L1 if Lp isequal to 1, the norm L2 if Lp is equal to 2 and so on, Lp being strictlyhigher than zero.
 6. The method as claimed in claim 5, wherein theanisotropy index A is calculated using the following relationship:$A = \sqrt{\sum\limits_{m = 1}^{M}\left( {\tau_{1,m}^{*2} + \tau_{2,m}^{*2}} \right)}$7. The method as claimed in claim 1, wherein: the vector u_(jk) is avector of

²\{(0,0)}, the modification T_(j,k) has: for d=2, the following matrixform: T j , k = k  [ cos  ( α j ) - sin  ( α j ) sin  ( α j ) cos ( α j ) ] and, for d=3, one of the following matrix forms or acombination of these matrix forms: T j , k = k  [ cos  ( α j ) - sin ( α j ) 0 sin  ( α j ) cos  ( α j ) 0 0 0 y k ] T j , k = k  [ cos ( α j ) 0 - sin  ( α j ) 0 y k 0 sin  ( α j ) 0 cos  ( α j ) ] T j ,k = k  [ y k 0 0 0 cos  ( α j ) - sin  ( α j ) 0 sin  ( α j ) cos ( α j ) ] the converted image being obtained by multiplying thecoordinates of the position of each pixel by the matrix T_(j,k), and theK-increments are calculated using only those pixels of the image thatoccupy a position m belonging to a set E, this set E containing onlythose positions m that already exist in the image I and that, whateverthe modification T_(j,k), after application of this modificationT_(j,k), occupy a position that also already exists in the image I, andfor which the position “m−T_(j,k)·p” occupies a position that alsoalready exists in the image I; the converted image I_(j,k) obtained atthe end of the calculation of this K-increment only containing pixelsthe positions of which belonged to the set E.
 8. The method as claimedin claim 7, wherein the p-variations calculated in step c) are quadraticvariations calculated using the following formula:${W_{j,k}\lbrack m\rbrack} = {\frac{1}{n_{E}}{\sum\limits_{m \in E}\left( {V_{j,k}\lbrack m\rbrack} \right)^{q}}}$where q=2 and n_(E) is the number of positions that belong to the set E.9. A method for automatically classifying digital images depending onthe anisotropy of their texture, this method comprising acquiring aplurality of images each formed from a plurality of pixels; wherein themethod further comprises automatically calculating, for each of theacquired images, a respective anisotropy index by means of a method asclaimed in claim 1, and classifying the acquired digital images using anautomatic classifier, depending on the anisotropy index calculated foreach of said images.
 10. A non-transitory computer-readable data storagemedium having encoded thereon instructions for performing a method asclaimed in claim 1, when these instructions are executed by a computer.11. A computer for implementing claim 1, said computer being programmedto execute the following steps: a) acquiring a digital image formed frompixels, each pixel being associated with a light intensity and with aposition in space

^(d), where d is a natural integer higher than or equal to two; b)automatically converting) the acquired image in order to obtain aconverted image I_(j,k), the conversion comprising applying amodification T_(j,k) of the image that makes each pixel of the acquiredimage rotate from one position to the next by an angle α_(j) about apoint or an axis and that enlarges or shrinks the image by a factorγ_(k), where α_(j)=arg(u_(jk)) and γ_(k)=|u_(jk)|², u_(jk) being avector that completely characterizes the modification T_(j,k), theindices j and k uniquely identifying the angle α_(j) and the factorγ_(k), respectively, then, for each converted image, calculating aK-increment V_(j,k)[m] for each pixel of position m of the convertedimage, this K-increment being calculated by applying a convolutionkernel v by means of the following formula:${V_{j,k}\lbrack m\rbrack} = {\sum\limits_{p \in {\lbrack{0,L}\rbrack}^{d}}{{v\lbrack p\rbrack} \cdot {Z\left\lbrack {m - {T_{j,k} \cdot p}} \right\rbrack}}}$wherein the product T_(j,k)·p corresponds to the application of themodification T_(j,k) to the pixel that initially had the position p inthe image I; the convolution kernel v performs linear filtering andpossesses a characteristic polynomial Q_(v)(z) and a finite support[0,L]^(d), v[p] being the value of the convolution kernel v for theposition p, the characteristic polynomial Q_(v)(z) being defined by thefollowing formula: ∀ z ∈ d , Q v  ( z ) = ∑ p ∈ [ 0 , L ] d  v  [ p ] z p and meeting the following condition: ∀a∈[0,K]^(d) such that |a|≤Kthen${\frac{\partial^{a}Q_{v}}{{\partial z_{1}^{a^{1}}}\mspace{11mu} \ldots \mspace{14mu} {\partial z_{d}^{a^{d}}}}\left( {1,\ldots \mspace{11mu},1} \right)} = 0$wherein L is an acquired vector of [0, N]^(d) that parametrizes thekernel v, N is a vector belonging to

^(d) that codes the size of the image and the components of which arestrictly positive natural integers; the constant K is a nonzero acquirednatural integer; z is a vector of components z₁, z₂, . . . , z_(d);z^(p) designates the monomial z₁ ^(p1)*z₂ ^(p2)* . . . *z_(d) ^(pd);∂^(|a|)Q_(v)/∂z₁ ^(a1) . . . ∂z_(d) ^(ad) is the partial derivative ofthe polynomial Q_(v)(z) with respect to the components of the vector z,the symbol ∂z_(i) ^(ai) indicating a differentiation of the polynomialQ_(v)(z) of order a_(i) with respect to the variable z_(i), where z_(i)is the i-th component of the vector z and a_(i) is the i-th component ofthe vector a, i being an integer index higher than or equal to 0 andlower than or equal to d; the step b) being executed with n_(j)different angles α_(j) and, for each angle α_(j), with at least twodifferent factors γ_(k), n_(j) being an integer higher than or equal to2 so as to obtain at least four different converted images I_(j,k); c)for each different converted image I_(j,k), calculating a p-variationW_(j,k) of this converted image from said calculated K-increments; d)estimating the terms β_(j) of the following statistical regression:log(|W_(j,k)|)=log(|u_(jk)|²)*H+β_(j)+ε_(j,k), wherein H is the Hurstexponent of the acquired image; ε_(j,k) is an error term of theregression the statistical properties of which are predetermined;wherein the computer is also programmed to execute the following steps:e) estimating the scalar coefficients τ_(m) of an even function τ(θ)defined in [0; 2π] that minimizes the following criterion C:$C = {\sum\limits_{j = 1}^{n_{j}}\left( {\beta_{j} - {\tau \star {\Gamma \left( \alpha_{j} \right)}}} \right)^{2}}$wherein β_(j) are the terms estimated in step d), τ(θ) is the functiondefined by the following relationship for any angle θ belonging to [0;2π]:${\tau (\theta)} = {\sum\limits_{m = 0}^{M}{\tau_{m}{f_{m}(\theta)}}}$wherein M is a constant integer number higher than an acquired datum,τ_(m) are the scalar coefficients of the function τ(θ), f_(m)(θ) are thefunctions of a basis of π-periodic functions defined in the interval [0;2π], Γ(θ) is the function defined by the following relationship:${\Gamma (\theta)} = {\int\limits_{{\mathbb{R}} +}{{{\hat{v}({\rho\theta})}}^{2}\rho^{{{- 2}H} - 1}d\; \rho}}$wherein {circumflex over (V)} is the discrete Fourier transform of thekernel v, H is the Hurst exponent of the acquired image, ρ is thevariable of integration, the symbol “*” designates the circularconvolution of the functions τ(θ) and Γ(θ), f) then, calculating,depending on the estimate of the scalar coefficients τ_(m), ananisotropy index that characterizes the anisotropy of the image, thisindex varying monotonically as a function of the statistical dispersionof the values of the function τ(θ) for θ varying between 0 and π.